SUBROUTINE Advec3d_X(m1,m2,m3, ia,iz,ja,jz,q0,u,den0,den1,dt,dxx,dd0,qn,mynum,nvar)
!-------------------------
  ! This subroutine calculates change in mixing ratio (Q0) during time
  !  step DT due to advection along a grid IDIM in length. Mixing ratios
  !  from host code (C) are loaded into Q0 array, which is updated to QN.
  !  Velocities (U) and fluxes (FLUX) are specified at cell FACES, having
  !  dimensions 0:IDIM. U, Q0, QN, DXX and FLUX indices defined here:
  !  Densities at beg, end time (DEN0, DEN1) defined in HOST CODE  
  !
  ! I grid->   |  1  |  2  |  I-1  |   I  |..   ..|  IDIM  | <- host grid
  ! U-array-> u(0)  u(1)  u(2)   u(i-1)  u(i)           u(IDIM)
  ! C-array->  | C(1)| C(2)| C(I-1)| C(I) |..   ..| C(IDIM)| mixing ratio
  ! DXX-arry-> | Dx1 | Dx2 | DxI-1 | DxI  |..   ..| DxIDIM |
  ! Density->  | Dd1 | Dd2 | DdI-1 | DdI  |..   ..| DdIDIM |
  !                 Q0 defined along 0 - IDIM+1 cells:
  !    |       | QN  | QN  |  QN   |  QN  |       |   QN   |        |
  
  !    |   Q0--|-Q0--|-Q0--|--Q0 --|--Q0--|..   ..|-- Q0 --|--Q0    |
  !    |    0  | 1   |  2  | I-1   |  I   |       |  IDIM  | IDIM+1 |
  !   lower BC |             <---   Q0 grid   --->         | upper BC
  !           Boundary conditions are stored in Q0 cells 0 & IDIM+1 
  !
  !  Input to this subroutine, provided in common /sub/, and the calling
  !  arguments to this subroutine:
  !     IDIM - #of grid cells being updated
  !     Q0(0:IDIM+1)- Initial mixing ratio along 1-D array, with two
  !                 additional boundary value mixing ratios padded into the
  !                 0th and IDIM+1 cell locations
  !     U(0:IDIM)- velocities BETWEEN grid cells (at the "higher-I" edges of
  !                each grid cell in the array, units consistent with DX, DT
  !     DEN0(IDIM)- Initial fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in Calling code
  !     DEN1(IDIM)- Updated fluid density, which needs to be updated during
  !                 multi-dimensional calculations, as noted in calling code
  !     DT-         time step- units consistent with U
  !     DXX(IDIM)-  Grid cell length along advection direction, Units
  !                   consistent with DT and U
  !     DD0(0:IDIM)- Initial fluid density flowing BETWEEN each grid cell
  !                  (remains constant for all dimensions at the initial
  !                  fluid density of the 1st dimension of a 2-3 D calculation
  !               one can use UPSTREAM density here (DD0(I)= RHO0(I) if u>0
  !               or DD0(I)= RHO0(I+1) if u<0) where RHO0 is the initial
  !               fluid density at the beginning of the 1st dimensional
  !               advection step of a 2 or 3 D advection calculation done one
  !               step at a time 
  !
  !  Output of this subroutine is an updated mixing ratio array QN(IDIM)
  !

  use node_mod, only : nodei0,nodej0
  
  use rpara, only: nmachs
       
  use mem_grid, only:  ngrids, nnxp,nnyp

!LFR>   use ParLib50, only: &
!LFR>        parf_wait_nostatus, &
!LFR>        parf_get_noblock, &
!LFR>        parf_send_noblock

  IMPLICIT none  
 
  INCLUDE "i8.h"
 
  INTEGER,INTENT(IN)                        :: m1,m2,m3, ia,iz,ja,jz,mynum,nvar
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3) :: q0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: u
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: den0
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: den1
  REAL   ,INTENT(in)                        :: dt
  REAL   ,INTENT(in),DIMENSION(m2,m3)       :: dxx
  REAL   ,INTENT(in),DIMENSION(m1,m2,m3)    :: dd0
  REAL   ,INTENT(OUT),DIMENSION(m1,m2,m3)   :: qn

  REAL,DIMENSION(m1,m2,m3)    :: flux
  REAL,DIMENSION(m1,m2,m3)    :: vcmax
  REAL,DIMENSION(m1,m2,m3)    :: vcmin
  
  LOGICAL,DIMENSION(m1,m2,m3) :: imxmn
  
  REAL,PARAMETER                       :: zr0=0.0
  
  INTEGER :: idime
  INTEGER :: i,j,k
  REAL    :: cf,cf1,ck1,ck2,x1,x1n,cs

  integer, parameter :: debug =0, dump=0
  integer, parameter :: on = 1, up=1,down=0
  
  INTEGER :: ii,ji,ii0,ji0,ie,je,ie0,je0,ipos,iia,iiz
  
  integer :: nf,itag
  

  REAL, DIMENSION(2*m1,ja:jz) :: buffer
   
  nf=80+mynum
  idime = m2 ! x dir

  qn=q0

  imxmn=.false.
  flux=0.0
  ! imxmn(idime+1)=.false.
  
  !WRITE (nf,FMT='(A)') 'X Begin'; CALL flush(nf)
  
  ii=ia+nodei0(mynum,ngrid)
  ii0=ia+nodei0(mynum,ngrid)-1
  ie=iz+nodei0(mynum,ngrid)
  ie0=iz+nodei0(mynum,ngrid)+1
  !print *,nvar

  DO j=ja,jz

    DO k=2,m1-1 
       ! Update mixing ratios and limit Fluxes going UP where u>0
       !  First assume upstream flux at edge of domain
       IF(u(k,1,j)>=zr0) flux(k,1,j)= q0(k,1,j)*u(k,1,j)*dt*dd0(k,1,j)
    END DO

    ji=j+nodej0(mynum,ngrid)
    IF(north(ngrid,mynum,ji)/=0) THEN
        itag=columnTag(ngrid,mynum,north(ngrid,mynum,ji),ji,1)
        IF (dump==on) WRITE (nf,FMT='("Get - Np: ",I2.2," J: ",I3.3," North: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,north(ngrid,mynum,ji),itag;CALL flush(nf)
        CALL get(buffer(:,j),2*m1,north(ngrid,mynum,ji),itag)
     END IF
   END DO 
  
   DO j=ja,jz
  
    ji=j+nodej0(mynum,ngrid)
    IF(north(ngrid,mynum,ji)/=0) THEN
       itag=columnTag(ngrid,mynum,north(ngrid,mynum,ji),ji,1)
       IF (dump==on)  WRITE (nf,FMT='("Wait - Np: ",I2.2," J: ",I3.3," Tag: ",I5.5)') &
	                     mynum,j,itag;CALL flush(nf)
       CALL parf_wait_nostatus(handle(itag))
       !IF (dump==on) WRITE (nf,FMT='(A)') 'Unbuffering';CALL flush(nf)
       CALL unPackBuffer(m1,buffer(:,j),flux(:,1,j),imxmn(:,1,j))
       IF (dump==on) CALL dumpData(m1,nf,flux(:,1,j),imxmn(:,1,j),buffer(:,j))
    END IF
   ! Identify local max and min, specify mixing ratio limits at new time
  DO  i=2,idime-1 ! ia,iz-1 or 1,iz-1

     DO k=2,m1-1 
				  
        imxmn(k,i,j)=q0(k,i,j)>=max(q0(k,i-1,j),q0(k,i+1,j)) .or. & !=true if local
              q0(k,i,j)<=min(q0(k,i-1,j),q0(k,i+1,j))       !       extrema
        ck1= q0(k,i,j)
        ck2= q0(k,i,j)
        if(u(k,i,j  )< zr0) ck1= q0(k,i+1,j)
        if(u(k,i-1,j)>=zr0) ck2= q0(k,i-1,j)
        vcmax(k,i,j)= max( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        vcmin(k,i,j)= min( q0(k,i,j), ck1, ck2 )                      ! Eq-7
        !  VCMAX and VCMIN are the absolute physical limits to the
        !     mixing ratio at t+dt. If these limits are ever violated,
        !     non-monotonic (oscillatory) behavior in solution results
        !

        IF(u(k,i,j)<zr0) CYCLE
     
        IF(u(k,i-1,j)<zr0) THEN
           flux(k,i,j)= q0(k,i,j)*u(k,i,j)*dt*dd0(k,i,j)    !  outflow-only cell     
        ELSE                              !      use upstream
           x1= dt*u(k,i,j)/dxx(i,j)               ! Courant number
           x1n= (1.-x1)*(q0(k,i+1,j)-q0(k,i-1,j))/4.
           !
           ! First, estimate mixing ratio in outflowing fluid (Cf)
           cf= q0(k,i,j) + x1n                                       !Eq-4a
           !
           !   Check to see if there is a peak (min) upwind and/or
           !    downwind of cell face
           IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
           IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +       (1.75 -.45*x1)*x1n   !Eq-10a
           !        CF= Q0(k,i,j) + 5.*X1N   ! uncomment this line for "full sharp"
           !
           !   Limit Cf to be between mixing ratio on either side of edge
           !      where flux is being calculated
           cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i+1,j))  ), MAX(q0(k,i,j),q0(k,i+1,j)) )
           !- for debug purposes only
	   if(debug==on)        CF1= CF     ! This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
           !                       implemented this subroutine constant mixing
           !                       ratios should be maintained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
        
	   !
           !   Calculate mixing ratio at new time, but limit to physically
           !    reasonable values
           qn(k,i,j)= MAX(vcmin(k,i,j),MIN(vcmax(k,i,j),          &   !eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-x1*cf1*dd0(k,i,j)+flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))

	   ! for debug purposes only
           if(debug==on)  QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-X1*CF1*DD0(k,i,j)+FLUX(k,I-1,j)/DXX(I,J))/DEN1(k,i,j)
           !            This statement IGNORES monotonic limitations.
           !                       you should uncomment this line and run your
           !                       advection calculation wth constant initial
           !                       mixing ratios everywhere. If you have properly
	   !                       implemented this subroutine constant mixing
           !                       ratios should be maintiained
           ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
           !   print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
	   !       if(qn(i) < 99.999) pause

	   !
           !   Re-calculate OUTFLOWING flux before moving on to next cell
           !    Flux = CF1*X1*DD0 but it must be adjusted if a monotonic limit
           !    is encountered.
           flux(k,i,j)= dxx(i,j)*(q0(k,i,j)*den0(k,i,j) - qn(k,i,j)*den1(k,i,j)) + flux(k,i-1,j)
        END IF                                                  !Eq-9a
     END DO
     
  END DO
END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell 1 if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !        IF(U(IDIM-1)>=ZR0.AND.U(IDIM)>=ZR0)
  !     &  QN(1)=(Q0(1)*DEN0(1)-FLUX(1)/DXX(1)+FLUX(IDIM)/DXX(1))/DEN1(1)
  !      END IF 
  !
  ! Update mixing ratios and limit Fluxes going DOWN where u<0
  !  The logic of this loop through the grid line is identical
  !  to the "DO 10" Loop above, only you start at the highest I
  !  edge and work backwards to I=1
  !


  DO j=ja,jz

     ji=j+nodej0(mynum,ngrid)
     IF(south(ngrid,mynum,ji)/=0) THEN
         itag=columnTag(ngrid,mynum,south(ngrid,mynum,ji),ji,2)
         CALL packBuffer(m1,flux(:,idime-1,j),imxmn(:,idime-1,j),buffer(:,j))
         IF (dump==on) WRITE (nf,FMT='("Put - Np: ",I2.2," J: ",I3.3," South: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,south(ngrid,mynum,ji),itag;CALL flush(nf)	 
	 IF (dump==on) CALL dumpData(m1,nf,flux(:,idime-1,j),imxmn(:,idime-1,j),buffer(:,j))
        CALL send(buffer(:,j),2*m1,south(ngrid,mynum,ji),itag)
     END IF  
  END DO

   !IF (nvar==3) THEN
   !    WRITE (nf,FMT='(A,2(I3.3,1X))')  'flux(1:2,:,7) &  imxmn(1:2,:,7)  --- ', &
    !                                    1+nodei0(mynum,ngrid),m2+nodei0(mynum,ngrid)
    !   DO k=1,2
!	  WRITE (nf,FMT='(I3.3,1X,$)') k!
!	  DO i=1,m2
!	     WRITE (nf,FMT='(E15.6,1X,L1,1X,$)') flux(k,i,7),imxmn(k,i,7)
 !         END DO
!	  WRITE (nf,FMT='A') ''
!       END DO
!       CALL flush(nf)
!    END IF


 DO j=ja,jz

     ji=j+nodej0(mynum,ngrid)
     IF(south(ngrid,mynum,ji)/=0) THEN
        itag=columnTag(ngrid,mynum,south(ngrid,mynum,ji),ji,1)
        IF (dump==on) WRITE (nf,FMT='("Get - Np: ",I2.2," J: ",I3.3," South: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,south(ngrid,mynum,ji),itag;CALL flush(nf)
        CALL get(buffer(:,j),2*m1,south(ngrid,mynum,ji),itag)
     ELSE
        DO k=2,m1-1       
           IF(u(k,idime-1,j)<zr0) flux(k,idime-1,j)= &
    	       q0(k,idime,j)*u(k,idime-1,j)*dt*dd0(k,idime-1,j)
        END DO
     END IF
  END DO 


  DO j=ja,jz
    ji=j+nodej0(mynum,ngrid)
    IF(south(ngrid,mynum,ji)/=0) THEN
	itag=columnTag(ngrid,mynum,south(ngrid,mynum,ji),ji,1)
          IF (dump==on) WRITE (nf,FMT='("Wait - Np: ",I2.2," J: ",I3.3," Tag: ",I5.5)') &
	                     mynum,j,itag;CALL flush(nf)
	  CALL parf_wait_nostatus(handle(itag))
          !WRITE (nf,FMT='(A)') 'Unbuffering';CALL flush(nf)
          CALL unPackBuffer(m1,buffer(:,j),flux(:,idime-1,j),imxmn(:,idime,j))
	  IF (dump==on) CALL dumpData(m1,nf,flux(:,idime-1,j),imxmn(:,idime,j),buffer(:,j))
     END IF
   !LFR-> from iz to ia


  !cycle!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<



   DO i=idime-1,2,-1 !iz,ia,-1
     ! print*,'part 2'
     DO k=2,m1-1 
        
        IF(u(k,i-1,j)>=zr0) THEN           ! Inflow-only cell
        
	   IF(u(k,i,j)<zr0) qn(k,i,j)=  MAX(  vcmin(k,i,j),   MIN(   vcmax(k,i,j),&
                  (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j) + flux(k,i-1,j)/dxx(i,j))/den1(k,i,j) ))
	ELSE
              x1=  dt*ABS(u(k,i-1,j))/dxx(i,j)     ! Courant number
              x1n= (1.-x1)*(q0(k,i-1,j)-q0(k,i+1,j))/4.
              cf= q0(k,i,j) + x1n                                       !Eq-4b
              IF(imxmn(k,i+1,j)) cf= q0(k,i,j) +MAX(1.5,1.2  +.6 *x1)*x1n   !Eq-10b
              IF(imxmn(k,i-1,j)) cf= q0(k,i,j) +   (1.75 -.45*x1)*x1n       !Eq-10a
              cf1= MIN( MAX( cf, MIN(q0(k,i,j),q0(k,i-1,j)) ), MAX(q0(k,i,j),q0(k,i-1,j)) )
              !- for debug purposes only
	      if(debug==on)  CF1= CF	! This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
              !			   ratios should be maintained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              
	      IF(u(k,i,j)>=zr0) cf1= q0(k,i,j)     ! outflow-only cell upstream
              qn(k,i,j)= MAX(  vcmin(k,i,j),  MIN(   vcmax(k,i,j), 	  &   !Eq-3&8
                    (q0(k,i,j)*den0(k,i,j)-flux(k,i,j)/dxx(i,j)-x1*cf1*dd0(k,i-1,j))/den1(k,i,j) ))
              !print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
          
	      !- for debug purposes only
              if(debug==on) QN(k,i,j)=(Q0(k,i,j)*DEN0(k,i,j)-FLUX(k,I,j)/DXX(I,J)-X1*CF1*DD0(k,I-1,j))/DEN1(k,i,j)
              !		This statement IGNORES monotonic limitations.
              !			   you should uncomment this line and run your
              !			   advection calculation wth constant initial
              !			   mixing ratios everywhere. If you have properly
              !			   implemented this subroutine constant mixing
	      !			   ratios should be maintiained
              ! DEBUG    DEBUG   DEBUG    DEBUG   DEBUG
              !  print*,'i q= ', i,qn(i),line,q0(i),q0(i-1),q0(i+1)
              !  print*,FLUX(I),x1,cf1!,dd0(i-1),den1(i)
              !  print*,q0(i)-FLUX(I)-x1*cf1!,dd0(i-1),den1(i)
	      !      if(qn(i) < 99.999) pause

	      flux(k,i-1,j)=dxx(i,j)*(qn(k,i,j)*den1(k,i,j) - q0(k,i,j)*den0(k,i,j)) + flux(k,i,j)!Eq-9b
           
	     !write(mynum,*) i+nodei0(mynum,ngrid),qn(k,i,j),q0(k,i,j),flux(k,i,j),flux(k,i-1,j); call flush(mynum)
             !if(i+nodei0(mynum,ngrid) == 16) then
   	     !  write(mynum,*)  vcmin(k,i,j), vcmax(k,i,j),q0(k,i,j),cf
   	     ! 	write(mynum,*)  x1,cf1,dd0(k,i-1,j),den0(k,i,j),den1(k,i,j)
             ! endif     ; call flush(mynum)
        END IF                                                  
     END DO

  END DO

  !END DO
  !
  ! If periodic boundary conditions are assumed, it is necessary
  !   to recalculate the updated mixing ratio at cell IDIM if there
  !   is inflow to that cell from the boundary between IDIM and 1
  !   Here these statements are commented out, but should be uncommented
  !   if this subroutine is needed for periodic boundary conditions,
  !   and then one of the calling arguements to the subroutine is IPERIOD
  !   which is set to "1" if you assume period boundary conditions
  !      IF(IPERIOD==1) THEN
  !      IF(U(1).LT.ZR0.AND.U(IDIM).LT.ZR0)
  !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
  !     &  QN(IDIM)=(Q0(IDIM)*DEN0(IDIM)-FLUX(0)/DXX(IDIM)+FLUX(IDIM-1)/
  !     &                DXX(IDIM))/DEN1(IDIM)
  !      END IF
  !
ENDDO !- big loop y-z

1000 continue


  DO j=ja,jz

    je=j+nodej0(mynum,ngrid)
    ji=j+nodej0(mynum,ngrid)

     !sending last column  (i=idime-1)  to east node
     IF(north(ngrid,mynum,ji)/=0) THEN
        itag=columnTag(ngrid,mynum,north(ngrid,mynum,ji),ji,2)
        CALL packBuffer(m1,flux(:,1,j),imxmn(:,2,j),buffer(:,j))
	IF (dump==on) WRITE (nf,FMT='("Put - Np: ",I2.2," J: ",I3.3," North: ",I2.2," Tag: ",I5.5)') &
	                     mynum,j,North(ngrid,mynum,ji),itag;CALL flush(nf)	
        CALL send(buffer(:,j),2*m1,north(ngrid,mynum,ji),itag)
	IF (dump==on) CALL dumpData(m1,nf,flux(:,1,j),imxmn(:,2,j),buffer(:,j))
     END IF  
  END DO
  !WRITE (nf,FMT='(A)') 'X ended'; CALL flush(nf)
  !STOP 500
END SUBROUTINE Advec3d_X
!--------------------------------------------------------------------------